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Via computer simulations of the standard binary Lennard- Jones glass former 
we have obtained in a systematic way a large set of close-by pairs of min- 
ima on the potential energy landscape, i.e. double-well potentials (DWP). 
We analyze this set of DWP in two directions. At low temperatures the 
symmetric DWP give rise to tunneling systems. We compare the resulting 
low-temperature anomalies with those, predicted by the standard tunneling 
model. Deviations can be traced back to the energy dependence of the relevant 
quantities like the number of tunneling systems. Furthermore we analyze the 
local structure around a DWP as well as the translational pattern during 
the transition between both minima. Local density anomalies are crucial for 
the formation of a tunneling system. Two very different kinds of tunneling 
systems are observed, depending on the type of atom (small or large) which 
forms the center of the tunneling system. In the first case the tunneling sys- 
tem can be interpreted as a single-particle motion, in the second case it is 
more collective. 

PACS numbers: 6143Fs 61.20Ja 



1. INTRODUCTION 

Almost all kinds of disordered solids show anomalous physical behavior 
at temperatures around 2 K Many of the observed features can be 
explained by the Standard Tunneling Model (STM)^Eland its generalization, 
the Soft-Potential Model EEEEI The basic idea of the STM is to postulate 
the existence of a broad distribution of Two Level States (TLS). A TLS 
can be represented by a single degree of freedom, moving in a double well 
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shaped potential, thereby acting as a bistable mode. The number of atoms, 
participating in this motion, does not enter in the model. According to the 
STM the distribution of TLS is chosen as 

P(A,A ) = P /A , (1) 

with an energy splitting 

E = v/A 2 + A§, (2) 

where A is the asymmetry and Ao the tunneling matrix element. The inde- 
pendence of Po from E is an important ingredient for the analytical predic- 
tions of the STM. The TLS couple to strain and electric fields and therefore 
influence the heat capacity, thermal conductivity, sound absorption, dielec- 
tric response and other quantities; see the review by Phillips^. 

The STM predicts a linear dependence of the heat capacity on tempera- 
ture and a quadratic dependence of the thermal conductivity on temperature. 
The observed small deviations from the predicted thermal conductivity be- 
havior can be explained by an energy dependence of Po an d the deformation 
potential, describing the coupling to strain Apart from these minor de- 
viations the STM gives a good general agreement with experimental results 
down to temperatures around 100 mK. 

Below this temperature, however, deviations from the STM occur^ 12 * 1 ^ . 
It is believed that this behavior is caused by interacting TLS which are not 
covered by the STM. The strong change of the dielectric response on an 
applied weak magnetic field which has been observed l 14 l 15 l 16 l are also not 
consistent with the STM. As shown recently, these effects can be under- 
stood by taking into account the interaction of the quadrupole moments of 
the nuclei, involved in the dynamics, with the local electric field gradient 
tensors 

Whereas it is very difficult to obtain microscopic information about the 
nature of TLS from experiments, such information can be supplied by com- 
puter simulations l 18 l 19 l 2 ° l 21 l. Most bistable modes will display an asymmetry 
much larger than 1 K. In general, energy differences between adjacent energy 
minima in glass-forming systems may range up to values much larger than 
T g Therefore, modes with asymmetries less than T g may be regarded as 
relatively symmetric, even if they do not contribute to the low-temperature 
anomalies. We will denote them as double-well potentials DWP. In this 
manuscript we perform a systematic search for DWP and analyze their mi- 
croscopic properties. This work extends previous work by Heuer and Silbey 
l 2 ^ 1 ^ in which circa 300 DWP could be identified for a binary Lennard- Jones 
system. This relatively small number of DWP was already sufficient to 
obtain direct information about the participation ratios and the absolute 
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number of DWP. In particular it was possible to extract the properties of 
the DWP with asymmetries less than 1 K, i.e. the TLS. In this way the 
deviations from the STM, the dependence of P on energy and its conse- 
quences could be elucidated. Due to advances in computer technology and 
application of a new algorithm to locate DWP we are now able to obtain a 
much larger set of DWP. In a recent paper we have shown that this set of 
DWP is not hampered by computational limitations like very fast cooling 
schedules as compared to the experiment The goal of the present work 
is twofold. First, we repeat our previous analysis to check the predictions of 
the STM. Second, we study the microscopic nature of the DWP in detail. 
This analysis requires a large set of DWP and was not possible before. 



2. COMPUTATIONAL DETAILS 

As a model system we have used a binary mixture Lennard-Jones sys- 
tem with 80% A-particles and 20% B-particles (BMLJ) V*WWV>W4 BMLJ 
is one of the standard glass-forming systems with very good glass-forming 
properties. The used potential is of the type 

V a/3 = A- e a(3 [(a af3 /r) 12 - (o a p/rf] + (a + b • r), (3) 

with cjab = O.80-AA, obb = O.880-AA, cab = l-5e A A, £bb = 0.5e A A, rn B = 
0.5mA- The simulation cell is a cube with a fixed edge length according to 
the number of particles and an exact particle density of D = 1.2. Periodic 
boundary conditions are used to minimize finite size effects and a linear 
function a + b ■ r has been added to the potential to ensure continuous 
energies and forces at the cutoff r c = 1.8. The units of length, mass and 
energy are oaa, mA, £AA> the time step within these units is set to 0.01. 

Previously, this potential has been used to mimic NiP E3 with 62 Ni and 
31 P, using o"aa = 2.2 A, the average mass per particle as 55.8 g/mol and 
€aa = 7765 J/mol. With this choice T=l in LJ-units corresponds to 934 
K. We just mention in passing that our density is the standard density for 
the BMLJ system, but is circa 20% higher than the density used for the 
mapping on the NiP system^. 

We apply molecular dynamics simulations using the velocity Verlet algo- 
rithm to equilibrate and generate a set of independent configurations. After 
minimizing these configurations we are systematically looking for nearby lo- 
cal energy minima and finally check whether these minima are connected by 
a simple barrier of first order with the respective starting minimum. Details 
of this search algorithm can be found in^J We have studied different sys- 
tem sizes (N = 65, 130, 195, 260) to identify possible finite size effects. For 
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the DWP properties, reported m*=Q, as well as for the results in the present 
work no finite size effects are present. The observation, that the structure 
of the DWP is not size dependent, even for our very small systems, is in 
accordance with data for the incoherent scattering function and the radial 
distribution function which already for N = 60 are close to the macroscopic 
limit 123. Furthermore it is known from previous work E3E2I that the number 
of participating particles is much smaller than the actual system size for 
the large majority of observed DWP. In what follows we restrict ourselves 
to system sizes N = 65 and N = 130. The smaller system sizes have the 
advantage that the efficiency of the DWP location algorithm is best ^ and 
to a good approximation a complete set of DWP is found, which enables us 
to obtain an estimate for the density of DWP. 

Of particular importance is the determination of the saddles between 
adjacent minima. The energy of the saddle is one of the major ingredients 
to calculate the tunneling matrix element. We have employed an algorithm, 
recently developed for the analysis of supercooled liquids above the glass 
transition ^3. 

Throughout this section we will use the following definitions for dis- 
tances. 



N 

d mw {fi,r 2 ) = J2 K> + d ly + d tz) ■ ^ ( 4 ) 



dmwrp — d mw 

! f 'trans. state ) + d mw (f t 

rans.statei fl ) (5) 



d%nw(f\,r2) is the mass weighted distance between two configurations, and 
dmwrp is the mass weighted reactions path approximation between two min- 
ima. The mass weighting has been introduced because the tunneling matrix 
element Ao depends on the mass weighted distance. The generated sets of 
DWP have been truncated in their parameters to guarantee the completeness 
of the search^. The maximum distance d mv/Tp is 0.8 and the asymmetry is 
limited to 0.5 if not stated otherwise. Within this parameter range we have 
found 6522 DWP from 10009 starting minima for the 65 particle system and 
2911 DWP from 3100 starting minima for N = 130, which is lowered by 
an inefficiency of the location algorithm for larger systems If not stated 
otherwise all presented data are generated from systems fully equilibrated at 
Tequil = 0.5, which is slightly above the critical mode-coupling temperature 

E2lr r = o.45. 
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Figure 1. Correlation of the averaged barrier height V and asymmetry A to 
the distance of the minima. Extra data are presented for comparison of the 
original data to those generated according to the W{ distribution. 

3. RESULTS 
3.1. Comparison with Experimental Findings 

Our first goal is to calculate the density of TLS as expressed by 

P eS (E) = J* dA J* dA ^5(E - V / A|TA^)P(A, A ). (6) 
Within the STM one has 

P eS (E) = P . (7) 

The value of P e g can be extracted from experimental data^. 

Every DWP is characterized by the triplet (A, d mwrp , V) where V is the 
barrier height. There exist major statistical correlations between all three 
quantities; see Fig. ^ For example DWP with small distances between the 
minima naturally display small asymmetries and barrier heights. Therefore 
direct extrapolation to nearly symmetric DWP, i.e. A < 1 K, is not possi- 
ble. To solve this problem we map every triplet (A,d mwrp , V) on a triplet 
(^2,^3,^4) such that the DWP, described by the fourth-order polynomial 



Epatix) = W2 * x 2 — W3 * x 3 + U>4 * x 4 , 



(8) 
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Figure 2. The Pi{wi) distributions as obtained from the DWP for N = 65. 

is characteriz ed by the same triplet (A, d mwrp , V). This procedure has been 
introduced in USED Thus our set of DWP is formally mapped on a distribu- 
tion P(w2, W3, W4). For this mapping each minimum of the DWP is taken 
as x = 0, so that two sets of wi result from one set of minima. 

Formally, the Wi can be viewed as the Taylor-expansion coefficients 
around the minima. In a disordered system the individual terms are a sum of 
very different contributions from the different pair interactions of the B JLM 
potential. Thus one may speculate that the Wi are statistically uncorrelated, 
i.e. 

-P(w 2 , W3, ^4) = p 2 (w 2 ) * P3(w 3 ) * £4(104). (9) 

Via a least square fit in the range of parameters, in which a systematic 
location of DWP has been performed, the Pi(wi) have been determined. 
They are shown in Fig. |2 Via a statistical procedure, as outlined in E2I we 
have checked that the Wi are indeed uncorrelated. 

Please note that the distribution Pi{wi) is very different to the density 
qi{wi) which denotes the number of DWP with a given value of Wi. The 
reason is that the sub- volume of the (u>2, W3, tu^-space which contains the 
DWP in the given parameter range, is highly non-trivial. Only in case that 
this sub-volume were a simple cube both distributions should be identical. 

Based on the distributions Pi(wi) it is possible to generate DWP with 
the correct statistics. Thus the generated DWP should have the same prop- 
erties as the original DWP. This is exemplified in Fig. ^ where we show the 
correlations of the DWP three parameters for the original set of DWP as 
well as for the generated set of DWP. No relevant deviations are present. 
This underlines the validity of our assumption of uncorrelated w^. 

Our main goal is to generate DWP with very small asymmetry in order 
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Figure 3. The distributions qi(wi) for A < 2K as calculated from the ob- 
served DWP. 

to obtain information about typical TLS. For perfectly symmetric DWP one 
has 

w 2 ■ w A /wl = 1/4. (10) 

However this relation alone does not tell which range in the W2, w 3, W4-space 
is relevant for nearly symmetric DWP. This question is crucial, since from 
Fig. [21 it is obvious that the p% and are ill-defined for small W3 and W4, 
respectively. Therefore in a first step we have generated DWP with A < 2K. 
For this subset we have calculated the (^-distributions. They are shown in 
Fig. El On their basis it is possible to estimate which parameter range 
is essential for the generation of nearly symmetric DWP. Fortunately, it 
turns out that the range of small 1113 and W4 is irrelevant. Thus it is indeed 
possible to use our original set of DWP to generate symmetric DWP, i.e. 
TLS with the correct statistics. This shows that the range of asymmetric 
DWP contains sufficient information about nearly symmetric DWP. This 
justified our method to use change the parameterization of the DWP by 
using the u>, rather than the original parameters (A, d mwrp , V). 

For the TLS, generated in this way, we estimate Ao with help of the 
Wentzel Kramer Brillouin (WKB) approximation. The validity of the one 
dimensional WKB approximation is discussed in for the case of TLS in 
simulated argon clusters. They observe only minor changes in the eigenfre- 
quencies during the transition, which makes the WKB approximation usable. 
We make use of the WKB approximation in a slightly different form than 
used by Phillips in^ 



An = E Q ■ e -^y/2mV'/n\ 



(11) 
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Figure 4. Energy dependency of P e ff and -Ptot- The data show the weighted 
and unweighted cumulative number of TLS per energy and volume. Both 
curves can roughly be approximated by a power law. 




Figure 5. The effective distance decreases with higher energy splitting, 
because large distances correlate with a high asymmetry and low Ao- The 
curve can roughly be approximated by a power law. 
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Instead of approximating the barrier shape by a rectangle we used an in- 
verted parabola. This gives the additional factor 7r/4 in the exponent. 
The primes indicate that the values have been changed by taking the har- 
monic ground level energies ^h(2w2/m) ' 5 of the relevant configurations, 
into account, i.e. minima and saddle. This lowers the barrier height V and 
decreases the distance d. By using the values for Ao, E and P(Ao,A) 
and according to Eqn. 13 with E/ks limited to 1 K we obtain P e ff = 
1.3 ± 0.5- 10 J m . This value is in good agreement with previous results 
^and the experimentally observed values for NiP 34 (2 • 10 46 J _1 m -3 ). This 
last value, however, might be biased by the electron contribution in NiP and 
is somewhat higher than commonly found values for other glass-forming sys- 
tems E31 ran gi n g from 0.5 • 10 45 J _1 m -3 to 3 • 10 45 J~ :L m~ 3 . In disagreement 
with the STM we observe an energy dependence of P e g(E) which can be 
characterized as P e g(E) oc E{ with 5\ = 0.14, see Fig. 4. 

In the next step we compare the temperature dependence of the ther- 
mal conductivity and the heat capacity calculated from simulated TLS with 
experimentally found values; see Heuer and Silbey ^ for the theoretical 
background. The thermal conductivity is experimentally found to behave 
like 

K (T)oc-^ocT 2 ^ (12) 

with (3 small and positive. From this and 

Pes oc E Sl (13) 

7 2 oc 4j(E) oc E*» (14) 
follows (3 = 5\ + 82 ■ d e g is defined as 



f > , , , , Ylid 2 • Af /E? 

d '" {E) = SAy E f • (15) 



where the primes indicate that the sums include all TLS i with an energy 
splitting smaller than E. The energy dependence of d e ff is shown in Fig. |S1 
yielding 82 = —0.17 and (3 ~ 0.0. For the calculation of the heat capacity 
experimentally limited relaxation times have to be taken into account. The 
heat capacity is experimentally found to be 

C(T) oc T 1,1-1 - 3 . (16) 

For a theoretical analysis one may use the expression 



2 



C(T,r cxp )oc r Pt ° t( f' T : x 2 p)E sech 
Jo 4k B T 



E 



2K B T 



dE. (17) 
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Figure 6. Dependence of the P e g on the Franck-Condon factor. N=65, 



Here Ptot is the total number of TLS per energy and volume with energy 
splitting below E and with relaxation times shorter than r exp . For a fixed 
relaxation time r exp = 5s we obtain P tot (E, r exp ) oc E 1+5s with 63 = 0.21; 
see Fig. 4. The value of ^3 is not very sensitive to the exact choice of r exp . In- 
serting this relation into Eg 1171 one obtains C(T,r exp ) oc T 1+<5a in agreement 
with the experiment. 

The comparison between experimentally determined results and simu- 
lated results has to be done with care as the tunneling matrix element is 
experimentally decreased by a mismatch of the vibrational modes of the two 
states giving rise to an additional factor 1251 



F is the Franck-Condon factor and does mainly depend on the spectral den- 
sity of eigenfrequencies and their coupling to the system. The simulated 
number of TLS is therefore not directly comparable with experimentally de- 
termined values. To check the possible influence of e~ F on our results some 
values, we have calculated P e g for different values of this factor; see Fig. |f)J 
Interestingly, P e g is slightly increased for 0.01 < e _F < 1. The increase is 
caused by the fact that TLS with E > IK are lowered in their splitting en- 
ergy such that the number of TLS below 1 K is increased. For even smaller 
values of the Franck-Condon factor, however, the Ag-factor in EqEl results 
in a significant reduction of P e d . 



E&1K 



exp — 



A 



(18) 
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Figure 7. Cosine of the angle between the transition vector from the first 
minimum configuration to the saddle configuration and the vector from the 
saddle to the next minimum configuration. The particles are sorted accord- 
ing to their total motion (i=0: fastest particle). It is distinguished whether 
central particle is a A- (triangles) or a B-particle (spheres). 




Figure 8. Sketch of the definition of the local environment of the central 
particle. 
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Figure 9. Left: Averaged A-particle distribution around the central B parti- 
cle. The filled circles are generated from a 65 particle system and the open 
circles are generated from a 130 particle system. There is no visible finite 
size effect. Right: Same as before but with a randomly chosen transition 
vector for the central particle. The density distribution corresponds to the 
radial distribution function. 
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Figure 11. G(r) of the minimized structures in comparison to special G(r), 
where only the central particles in interaction with A-particles are consid- 
ered. 



3.2. Microscopic Nature of the Two Level States 

In our recent publication ^ we have analyzed the particle which shows 
the largest displacement during the transition between both minima. It is 
denoted central particle. DWP with an A-particle as a central particle are 
denoted A-type DWP. B-type DWP are defined in analogy. Surprisingly, it 
turned out that 90% of all DWP are B-type DWP although only 20% of 
all particles are B-particle. Furthermore we see that for B-type DWP the 
displacement of the central particles was significantly larger than that of 
all other particles. For A-type DWP the distribution of displacements was 
more continuous. Stated differently, the transition in A-type DWP is more 
collective. 

In this subsection we elucidate the dynamics of DWP much closer. 
First we analyze whether the displacement between both minima is along a 
straight line in configuration space or whether it is curved. More specifically 
we calculate the angle between the transition vector from the first minimum 
to the saddle and from the saddle to the second minimum; see Fig. The 
results are sorted with the respect to the displacement of the different parti- 
cles (i=0: particle with the largest displacement, i.e. central particle; i=64: 
particle with the smallest displacement). For B-type DWP the central par- 
ticle basically moves along a straight path whereas the other particles move 
along a strongly curved path. Thus the DWP can be basically characterized 
as a one-particle motion of a B-particle. The other particles just support 
this transition by some complicated curved trajectory to optimize the en- 
ergy. The behavior is very different for A-type DWP where most particles 
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Figure 12. Particle distribution around the particle with the second-largest 
displacement (N=65). In analogy to Fig. Eland Fig. ^] the transition vector 
of this particle defines the y-axis. The filled circles mark the contribution 
from the central particle. 



move along a relatively straight line. Since the central particle is only little 
different in terms of its displacement (see above) it is not surprising that 
the central particle behaves similarly as compared to the other particles. 
Obviously, the collective dynamics in A-type DWP is realized by a rather 
straight translation of most particles. Note, however, that for both types 
of DWP there is a tendency to move in a less curved way for particles with 
larger displacements, even if the reaction path approximation is considered. 
A priori this is not a necessity. 

In the next step we take a closer look into the local distribution of par- 
ticles around the central particle. We therefore define the relative density 
in a cylinder around the central particle as shown in Fig. El The transition 
vector of the central particle defines the y-axis. Furthermore the distance 
vectors between the central particles and the other particles in the starting 
minimum are represented in the xy-plane. By taking into account appropri- 
ate phase space factors the situation of an ideal gas would correspond to a 
homogeneous density distribution in the xy-plane. Finally, we average over 
all DWP, distinguishing A-type and B-type DWP. We only consider the den- 
sity of A-particles around the central particles because the B-particles are 
of minor statistical relevance. Fig. |2l shows the results for the B-type DWP, 
which dominate the investigated BMLJ-system. The first striking observa- 
tion is the distinct structure (left), which is different from the average pair 
correlation function (right). The main structural feature is the very small 
density for (x w 0,y k <tab = 0.88). This directly shows that the central 
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Figure 13. Averaged projected displacement of all A-particles relative to the 
transition of the central particle. Results for N=65 and N=130 are shown. 



B-particles moves in a direction where free space is available. An increased 
density in the nearest neighbor shell (nn-shell) is found orthogonal to the 
direction of the translation vector. Thus the central particle jumps through 
this structure to get to the new minimum. Note that results for N=65 and 
N=130 are identical within numerical uncertainties. For the larger system 
one can see that also the second nearest neighbor shell reflects the proper- 
ties of nn-shell. On a qualitative basis the same properties are observed for 
A-type DWP; see Fig. I1UI It is, however, much less pronounced. This is a 
consequence of the fact that A-type DWP are more collective. 

A priori it is not clear whether this observation implies an anomalous 
structure around the central particle. This is checked by analyzing the radial 
distribution function G(r) for the average minimized structures and G(r) for 
the central particles of the minimized structures, see Fig. 1111 We can clearly 
see that for both types of DWP density is taken from the first nn-shell and 
transferred to a new peak between the first and the second nearest neighbor 
shell. Thus the observed structure in the local environment does not arise 
from the fact that the central transition vector chooses a special direction in 
an average environment, but that the environment is indeed altered. 

To illustrate the above observations we show a selected DWP which ex- 
plicitly displays the structural patterns; see Fig. 1141 Please note the particles 
in the intermediate shell, as well as the particles orthogonal to the transition 
vector. 

To go further into structural details we repeat the analysis, performed 
so far for the central particle, also for the particle with the second-largest 
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Figure 14. Type B-DWP as observed from the N=65 simulation. The small 
black central particle is a B-particle. The shown bonds have a length around 
0.9 the non bonded atoms have a distance to the B-particle between 1.11 
and 1.21. The two configurations correspond to the minima of a DWP. 



displacement. The data are presented in Fig.^Jfor both types of DWP. The 
structural changes for the B-type DWP are dramatic if compared to Fig. |H1 
For the B-type DWP the first shell is almost back to normal except for an 
increased probability to find particles just in front of the second particle. 
This density is mainly caused by the central particle. This means that the 
particle with the second-largest displacement mainly follows the central B- 
particle. The A-type DWP still shows a similar picture as Fig. llUf left) but 
much less pronounced, the central particle shows only a small tendency to 
occupy positions in front of the second particle (tiny filled circles ahead of 
the second particle). This again is a direct consequence of the fact that for 
A-type DWP the central particle is not behaving very different to the other 
participating particles. 

So far we have analyzed the structure around the central particle. In 
the final step we also elucidate the displacement of the particles around the 
central particle. The representation in Fig. EH is analog to those in Fig. |H1 
the x-y-plane is still the same, but now we do not analyze where particles are 
in the plane but how they move relative to the central transition vector. The 
displacement is projected into the x-y-plane. Both graphs in Fig. ^3 look 
similar, the differences between A- and B-type DWP are on a quantitative 
level. The general explanation for the observed displacements are simple. 
The fastest particle makes a jump from one minimum to the other and all 
other particles follow as if they were connected to the central particle by a 
spring, reflecting the attractive part of the Lennard-Jones potential. 
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4. SUMMARY 

We were able to present a huge set of DWP, which were found on a 
systematic basis. The number of TLS is consistent with experimental obser- 
vations and the deviations from the predictions of the STM can be explained 
in terms of the energy dependence of P e g and Ptot • These results were already 
part of older work ^ but are now reproduced on a much better statistical 
basis. 

With the analysis of the saddle angles during the transition and the 
structure of the second particle we underline the single particle nature of the 
transition in B-type DWP, as already proposed in^. As the B-type DWP are 
predominant in the investigated BMLJ-system it follows that single particle- 
type DWP determine the low temperature properties of BMLJ. The more 
collective A-type transitions are not as important for our system but can be 
used as a model for other substances, where all molecules have approximately 
the same size and mobility. 

We were also able to present a detailed analysis of the structure and 
the displacements of the DWP. It is shown that the presence of holes in the 
structure is a prerequisite for the formation of DWP. For B-type DWP one 
basically has a one-particle displacement. The other particles mainly follow 
the central particle. In contrast, for A-type DWP the dynamics are much 
more collective. 

In a next step we aim to compare the properties with that of network 
glasses like silicate in order to see how the structural differences are reflected 
in the low-temperature properties of glasses. 
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